# PoP - Policing Socio-Geographic Boundaries and Inequality
# script for creating Figures 2 & 3
# Figure 2 - Influence of Logged Crime on Logged Felony Arrests (Standardized), Conditional on Racial Boundary Status.
# Figure 3 - Influence of Logged Crime on Logged Misdemeanor Arrests (Standardized), Conditional on Racial Boundary Status.


# load packages
suppressPackageStartupMessages(
  
  {
    library(dplyr)
    library(tidyverse)
    library(ggplot2)
    library(haven)
    library(readxl)
    library(readr)
    library(areal)
    library(car)
    library(estimatr)
    library(magrittr)
    library(texreg)
    library(sandwich)
    library(jtools)
  }
)


# load and format all data
# First ATL
# Load Atlanta final dataset
load("atl_final.RData")

# Set up DVs
# log and +1 to variables (pop already logged)
atl_fin = atl_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property + 1),
         lviolentcrime = log(crime_violent + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
atl_fin$larrest_sd <- atl_fin$larrests - (mean(atl_fin$larrests)/sd(atl_fin$larrests))
atl_fin$lmisdemeanors_sd <- atl_fin$lmisdemeanors - (mean(atl_fin$lmisdemeanors)/sd(atl_fin$lmisdemeanors))
atl_fin$lfelonies_sd <- atl_fin$lfelonies - (mean(atl_fin$lfelonies)/sd(atl_fin$lfelonies))
atl_fin$lnonviolent_sd <- atl_fin$lnonviolent - (mean(atl_fin$lnonviolent)/sd(atl_fin$lnonviolent))
atl_fin$lviolent_sd <- atl_fin$lviolent - (mean(atl_fin$lviolent)/sd(atl_fin$lviolent))
atl_fin$lsociety_sd <- atl_fin$lsociety - (mean(atl_fin$lsociety)/sd(atl_fin$lsociety))

atl_fin$lperson_sd <- atl_fin$lperson - (mean(atl_fin$lperson)/sd(atl_fin$lperson))
atl_fin$lproperty_sd <- atl_fin$lproperty - (mean(atl_fin$lproperty)/sd(atl_fin$lproperty))

# now create binned racial blv measures
atl_fin <- atl_fin %>% mutate(boundary_quart = quantile(atl_fin$p_race_white_blv, prob=.75, na.rm=TRUE),
                              boundary_quart_dummy = ifelse(p_race_white_blv > boundary_quart,1,0 ))


# For Austin
# Load Austin final dataset
load("aus_final.RData")


# Set up DVs
# log and +1 to variables (pop already logged)
aus_fin = aus_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property + 1),
         lviolentcrime = log(crime_violent + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
aus_fin$larrest_sd <- aus_fin$larrests - (mean(aus_fin$larrests)/sd(aus_fin$larrests))
aus_fin$lmisdemeanors_sd <- aus_fin$lmisdemeanors - (mean(aus_fin$lmisdemeanors)/sd(aus_fin$lmisdemeanors))
aus_fin$lfelonies_sd <- aus_fin$lfelonies - (mean(aus_fin$lfelonies)/sd(aus_fin$lfelonies))
aus_fin$lnonviolent_sd <- aus_fin$lnonviolent - (mean(aus_fin$lnonviolent)/sd(aus_fin$lnonviolent))
aus_fin$lviolent_sd <- aus_fin$lviolent - (mean(aus_fin$lviolent)/sd(aus_fin$lviolent))
aus_fin$lsociety_sd <- aus_fin$lsociety - (mean(aus_fin$lsociety)/sd(aus_fin$lsociety))
aus_fin$lperson_sd <- aus_fin$lperson - (mean(aus_fin$lperson)/sd(aus_fin$lperson))
aus_fin$lproperty_sd <- aus_fin$lproperty - (mean(aus_fin$lproperty)/sd(aus_fin$lproperty))

# now create binned racial blv measures
aus_fin <- aus_fin %>% mutate(boundary_quart = quantile(aus_fin$p_race_white_blv, prob=.75, na.rm=TRUE),
                              boundary_quart_dummy = ifelse(p_race_white_blv > boundary_quart,1,0 ))


# For Boston
## Load Boston data
load("bos_final.RData")


# Set up DVs
# log and +1 to variables (pop already logged)
bos_fin = bos_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime + 1))



# create scaled DVs (to account for different pop. sizes and density in blocks)
bos_fin$larrests_sd <- bos_fin$larrests - (mean(bos_fin$larrests)/sd(bos_fin$larrests))
bos_fin$lmisdemeanors_sd <- bos_fin$lmisdemeanors - (mean(bos_fin$lmisdemeanors)/sd(bos_fin$lmisdemeanors))
bos_fin$lfelonies_sd <- bos_fin$lfelonies - (mean(bos_fin$lfelonies)/sd(bos_fin$lfelonies))
bos_fin$lnonviolent_sd <- bos_fin$lnonviolent - (mean(bos_fin$lnonviolent)/sd(bos_fin$lnonviolent))
bos_fin$lviolent_sd <- bos_fin$lviolent - (mean(bos_fin$lviolent)/sd(bos_fin$lviolent))
bos_fin$lsociety_sd <- bos_fin$lsociety - (mean(bos_fin$lsociety)/sd(bos_fin$lsociety))
bos_fin$lperson_sd <- bos_fin$lperson - (mean(bos_fin$lperson)/sd(bos_fin$lperson))
bos_fin$lproperty_sd <- bos_fin$lproperty - (mean(bos_fin$lproperty)/sd(bos_fin$lproperty))


# now create binned racial blv measures
bos_fin <- bos_fin %>% mutate(boundary_quart = quantile(bos_fin$p_race_white_blv, prob=.75, na.rm=TRUE),
                              boundary_quart_dummy = ifelse(p_race_white_blv > boundary_quart,1,0 ))




# For Chicago
## Load Chicago data
load("chi_final.RData")


# Set up DVs
# log and +1 to variables (pop already logged)
chi_fin = chi_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(property_crime + 1),
         lviolentcrime = log(violent_crime + 1))


# create scaled DVs (to account for different pop. sizes and density in blocks)
chi_fin$larrest_sd <- chi_fin$larrests - (mean(chi_fin$larrests)/sd(chi_fin$larrests))
chi_fin$lmisdemeanors_sd <- chi_fin$lmisdemeanors - (mean(chi_fin$lmisdemeanors)/sd(chi_fin$lmisdemeanors))
chi_fin$lfelonies_sd <- chi_fin$lfelonies - (mean(chi_fin$lfelonies)/sd(chi_fin$lfelonies))
chi_fin$lnonviolent_sd <- chi_fin$lnonviolent - (mean(chi_fin$lnonviolent)/sd(chi_fin$lnonviolent))
chi_fin$lviolent_sd <- chi_fin$lviolent - (mean(chi_fin$lviolent)/sd(chi_fin$lviolent))
chi_fin$lsociety_sd <- chi_fin$lsociety - (mean(chi_fin$lsociety)/sd(chi_fin$lsociety))
chi_fin$lperson_sd <- chi_fin$lperson - (mean(chi_fin$lperson)/sd(chi_fin$lperson))
chi_fin$lproperty_sd <- chi_fin$lproperty - (mean(chi_fin$lproperty)/sd(chi_fin$lproperty))

# now create binned racial blv measures
chi_fin <- chi_fin %>% mutate(boundary_quart = quantile(chi_fin$p_race_white_blv, prob=.75, na.rm=TRUE),
                              boundary_quart_dummy = ifelse(p_race_white_blv > boundary_quart,1,0 ))




# For Milwaukee
## Load Milwaukee data
load("mil_final.RData")


# Set up DVs
# log and +1 to variables (pop already logged) 
mil_fin = mil_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property + 1),
         lviolentcrime = log(crime_violent + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
mil_fin$larrest_sd <- mil_fin$larrests - (mean(mil_fin$larrests)/sd(mil_fin$larrests))
mil_fin$lmisdemeanors_sd <- mil_fin$lmisdemeanors - (mean(mil_fin$lmisdemeanors)/sd(mil_fin$lmisdemeanors))
mil_fin$lfelonies_sd <- mil_fin$lfelonies - (mean(mil_fin$lfelonies)/sd(mil_fin$lfelonies))
mil_fin$lnonviolent_sd <- mil_fin$lnonviolent - (mean(mil_fin$lnonviolent)/sd(mil_fin$lnonviolent))
mil_fin$lviolent_sd <- mil_fin$lviolent - (mean(mil_fin$lviolent)/sd(mil_fin$lviolent))
mil_fin$lsociety_sd <- mil_fin$lsociety - (mean(mil_fin$lsociety)/sd(mil_fin$lsociety))
mil_fin$lperson_sd <- mil_fin$lperson - (mean(mil_fin$lperson)/sd(mil_fin$lperson))
mil_fin$lproperty_sd <- mil_fin$lproperty - (mean(mil_fin$lproperty)/sd(mil_fin$lproperty))

# now create binned racial blv measures
mil_fin <- mil_fin %>% mutate(boundary_quart = quantile(mil_fin$p_race_white_blv, prob=.75, na.rm=TRUE),
                              boundary_quart_dummy = ifelse(p_race_white_blv > boundary_quart,1,0 ))


# For Seattle
## Load Seattle data
load("sea_final.RData")


# Set up DVs
# log and +1 to variables (pop already logged) 
sea_fin = sea_fin %>% 
  mutate(lmhhi = log(mhhi + 1),
         larrests = log(total_arrests + 1),
         lmisdemeanors = log(misdemeanor_arrests + 1),
         lfelonies = log(felony_arrests + 1),
         lnonviolent = log(nonviolent_arrests + 1),
         lviolent = log(violent_arrests + 1),
         lsociety = log(society_arrests + 1),
         lperson = log(person_arrests + 1),
         lproperty = log(property_arrests + 1),
         lcrime = log(crime_all + 1),
         lpropertycrime = log(crime_property + 1),
         lviolentcrime = log(crime_violent + 1))

# create scaled DVs (to account for different pop. sizes and density in blocks)
sea_fin$larrest_sd <- sea_fin$larrests - (mean(sea_fin$larrests)/sd(sea_fin$larrests))
sea_fin$lmisdemeanors_sd <- sea_fin$lmisdemeanors - (mean(sea_fin$lmisdemeanors)/sd(sea_fin$lmisdemeanors))
sea_fin$lfelonies_sd <- sea_fin$lfelonies - (mean(sea_fin$lfelonies)/sd(sea_fin$lfelonies))
sea_fin$lnonviolent_sd <- sea_fin$lnonviolent - (mean(sea_fin$lnonviolent)/sd(sea_fin$lnonviolent))
sea_fin$lviolent_sd <- sea_fin$lviolent - (mean(sea_fin$lviolent)/sd(sea_fin$lviolent))
sea_fin$lsociety_sd <- sea_fin$lsociety - (mean(sea_fin$lsociety)/sd(sea_fin$lsociety))
sea_fin$lperson_sd <- sea_fin$lperson - (mean(sea_fin$lperson)/sd(sea_fin$lperson))
sea_fin$lproperty_sd <- sea_fin$lproperty - (mean(sea_fin$lproperty)/sd(sea_fin$lproperty))

# now create binned racial blv measures
sea_fin <- sea_fin %>% mutate(boundary_quart = quantile(sea_fin$p_race_white_blv, prob=.75, na.rm=TRUE),
                              boundary_quart_dummy = ifelse(p_race_white_blv > boundary_quart,1,0 ))


###### moderation between race boundary dummy and total crime; on misdemeanor arrests ###


# atlanta, out of sample
mdata <- atl_fin %>% dplyr::select(lmisdemeanors_sd,boundary_quart_dummy,
                                   lcrime,ses_blv,lpop,p_race_white,age_15_35_male,
                                   hhi,lmhhi,pown,p_poverty,p_emp_unemployed,pcol,BG_CODE,pop,
                                   phys_bound, com_density) %>% na.omit()


atl1 = lm_robust(lmisdemeanors_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                   hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density, 
                 data = mdata, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

summary(mdata$lcrime)
low_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(0,length(lcrime)),
  lcrime = seq(0,6.9,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

library(prediction)
low_preds<- data.frame(cbind(prediction(atl1, low_dat, interval = "confidence", level=.95), low_dat))

head(low_preds)
hi_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(1,length(lcrime)),
  lcrime = seq(0,6.9,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

hi_preds<- data.frame(cbind(prediction(atl1, hi_dat, interval = "confidence", level=.95), hi_dat))

boundary<-rep("Boundary",length(low_preds))
ylow<-data.frame(low_preds)
ylow<-ylow %>% mutate(boundary="No Boundary")
yhi<-data.frame(hi_preds)
yhi<-yhi %>% mutate(boundary="Boundary")
preds<-data.frame(rbind(ylow,yhi))

atl_preds <- preds %>% mutate(city = "Atlanta")



# Austin, out of sample
mdata <- aus_fin %>% dplyr::select(lmisdemeanors_sd,boundary_quart_dummy,
                                   lcrime,ses_blv,lpop,p_race_white,age_15_35_male,
                                   hhi,lmhhi,pown,p_poverty,p_emp_unemployed,pcol,BG_CODE,pop,
                                   phys_bound, com_density) %>% na.omit()


aus1 = lm_robust(lmisdemeanors_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                   hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density, 
                 data = mdata, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

summary(mdata$lcrime)
low_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(0,length(lcrime)),
  lcrime = seq(0,8.1,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

low_preds<- data.frame(cbind(prediction(aus1, low_dat, interval = "confidence", level=.95), low_dat))

head(low_preds)
hi_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(1,length(lcrime)),
  lcrime = seq(0,8.1,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

hi_preds<- data.frame(cbind(prediction(aus1, hi_dat, interval = "confidence", level=.95), hi_dat))

boundary<-rep("Boundary",length(low_preds))
ylow<-data.frame(low_preds)
ylow<-ylow %>% mutate(boundary="No Boundary")
yhi<-data.frame(hi_preds)
yhi<-yhi %>% mutate(boundary="Boundary")
preds<-data.frame(rbind(ylow,yhi))

aus_preds <- preds %>% mutate(city = "Austin")


# boston, out of sample

mdata <- bos_fin %>% dplyr::select(lmisdemeanors_sd,boundary_quart_dummy,
                                   lcrime,ses_blv,lpop,p_race_white,age_15_35_male,
                                   hhi,lmhhi,pown,p_poverty,p_emp_unemployed,pcol,BG_CODE,pop,
                                   phys_bound, com_density) %>% na.omit()


bos1 = lm_robust(lmisdemeanors_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                   hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density, 
                 data = mdata, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

summary(mdata$lcrime)
low_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(0,length(lcrime)),
  lcrime = seq(0,7,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

low_preds<- data.frame(cbind(prediction(bos1, low_dat, interval = "confidence", level=.95), low_dat))

head(low_preds)
hi_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(1,length(lcrime)),
  lcrime = seq(0,7,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

hi_preds<- data.frame(cbind(prediction(bos1, hi_dat, interval = "confidence", level=.95), hi_dat))

boundary<-rep("Boundary",length(low_preds))
ylow<-data.frame(low_preds)
ylow<-ylow %>% mutate(boundary="No Boundary")
yhi<-data.frame(hi_preds)
yhi<-yhi %>% mutate(boundary="Boundary")
preds<-data.frame(rbind(ylow,yhi))

bos_preds <- preds %>% mutate(city = "Boston")


# chicago, out of sample
mdata <- chi_fin %>% dplyr::select(lmisdemeanors_sd,boundary_quart_dummy,
                                   lcrime,ses_blv,lpop,p_race_white,age_15_35_male,
                                   hhi,lmhhi,pown,p_poverty,p_emp_unemployed,pcol,BG_CODE,pop,
                                   phys_bound, com_density) %>% na.omit()


chi1 = lm_robust(lmisdemeanors_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                   hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + 
                   com_density, 
                 data = mdata, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

summary(mdata$lcrime)
low_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(0,length(lcrime)),
  lcrime = seq(0,9.25,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

low_preds<- data.frame(cbind(prediction(chi1, low_dat, interval = "confidence", level=.95), low_dat))

head(low_preds)
hi_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(1,length(lcrime)),
  lcrime = seq(0,9.25,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

hi_preds<- data.frame(cbind(prediction(chi1, hi_dat, interval = "confidence", level=.95), hi_dat))

boundary<-rep("Boundary",length(low_preds))
ylow<-data.frame(low_preds)
ylow<-ylow %>% mutate(boundary="No Boundary")
yhi<-data.frame(hi_preds)
yhi<-yhi %>% mutate(boundary="Boundary")
preds<-data.frame(rbind(ylow,yhi))

chi_preds <- preds %>% mutate(city = "Chicago")

## milwaukee out of sample

mdata <- mil_fin %>% dplyr::select(lmisdemeanors_sd,boundary_quart_dummy,
                                   lcrime,ses_blv,lpop,p_race_white,age_15_35_male,
                                   hhi,lmhhi,pown,p_poverty,p_emp_unemployed,pcol,BG_CODE,pop,
                                   phys_bound, com_density) %>% na.omit()


mil1 = lm_robust(lmisdemeanors_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                   hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound +
                   com_density, 
                 data = mdata, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

summary(mdata$lcrime)
low_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(0,length(lcrime)),
  lcrime = seq(0,7.5,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

low_preds<- data.frame(cbind(prediction(mil1, low_dat, interval = "confidence", level=.95), low_dat))

head(low_preds)
hi_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(1,length(lcrime)),
  lcrime = seq(0,7.5,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

hi_preds<- data.frame(cbind(prediction(mil1, hi_dat, interval = "confidence", level=.95), hi_dat))

boundary<-rep("Boundary",length(low_preds))
ylow<-data.frame(low_preds)
ylow<-ylow %>% mutate(boundary="No Boundary")
yhi<-data.frame(hi_preds)
yhi<-yhi %>% mutate(boundary="Boundary")
preds<-data.frame(rbind(ylow,yhi))

mil_preds <- preds %>% mutate(city = "Milwaukee")

## seattle, out of sample

mdata <- sea_fin %>% dplyr::select(lmisdemeanors_sd,boundary_quart_dummy,
                                   lcrime,ses_blv,lpop,p_race_white,age_15_35_male,
                                   hhi,lmhhi,pown,p_poverty,p_emp_unemployed,pcol,BG_CODE,pop,
                                   phys_bound, com_density) %>% na.omit()


sea1 = lm_robust(lmisdemeanors_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                   hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density, 
                 data = mdata, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

summary(mdata$lcrime)
low_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(0,length(lcrime)),
  lcrime = seq(0,7.4,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

low_preds<- data.frame(cbind(prediction(sea1, low_dat, interval = "confidence", level=.95), low_dat))

head(low_preds)
hi_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(1,length(lcrime)),
  lcrime = seq(0,7.4,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

hi_preds<- data.frame(cbind(prediction(sea1, hi_dat, interval = "confidence", level=.95), hi_dat))

boundary<-rep("Boundary",length(low_preds))
ylow<-data.frame(low_preds)
ylow<-ylow %>% mutate(boundary="No Boundary")
yhi<-data.frame(hi_preds)
yhi<-yhi %>% mutate(boundary="Boundary")
preds<-data.frame(rbind(ylow,yhi))

sea_preds <- preds %>% mutate(city = "Seattle")


#### stack 'em together and graph ####

# out of sample predictions graph

### out of sample predictions

atl_preds <- atl_preds %>% dplyr::select(lcrime,fitted.fit,
                                         fitted.lwr,fitted.upr,boundary,city)

aus_preds <- aus_preds %>% dplyr::select(lcrime,fitted.fit,
                                         fitted.lwr,fitted.upr,boundary,city)

bos_preds <- bos_preds %>% dplyr::select(lcrime,fitted.fit,
                                         fitted.lwr,fitted.upr,boundary,city)

chi_preds <- chi_preds %>% dplyr::select(lcrime,fitted.fit,
                                         fitted.lwr,fitted.upr,boundary,city)

mil_preds <- mil_preds %>% dplyr::select(lcrime,fitted.fit,
                                         fitted.lwr,fitted.upr,boundary,city)

# lou_preds <- lou_preds %>% dplyr::select(lcrime,fitted.fit,
                                         # fitted.lwr,fitted.upr,boundary,city)

sea_preds <- sea_preds %>% dplyr::select(lcrime,fitted.fit,
                                         fitted.lwr,fitted.upr,boundary,city)

d <- rbind(atl_preds, aus_preds, bos_preds, chi_preds,
           mil_preds,sea_preds)
head(d)

p <- ggplot(data=d, aes(x=lcrime, y=fitted.fit,group=boundary)) +
  facet_wrap(~ city, nrow = 2, scales = "free")+
  #geom_point(aes(shape = Boundary), color = "grey",alpha=.5)+
  geom_line(aes(y=fitted.fit, linetype=boundary, color = boundary))+
  geom_ribbon(aes(ymin=fitted.lwr, ymax=fitted.upr,linetype=boundary, color=boundary),alpha=0.1)+
  labs(y="Misdemeanor Arrests (Logged)",x="Crime (Logged)",
       title="Police Attention to Crime in and out of Boundary Zones")+
  theme_bw(base_size=12,base_family="Times")+
  theme(panel.border=element_rect(fill=NA, colour=NA), 
        legend.title = element_blank(),
        legend.position="bottom",
        panel.grid.minor=element_line(colour=NA))+
  theme(strip.background = element_rect(fill="white"))

ggsave(plot = p, filename = "figure2_plots.png",
       height = 8, width = 8)

##############
## Felonies ##
##############
names(atl_fin)
## Atlanta

# expected values #
mdata <- atl_fin %>% dplyr::select(lfelonies_sd,boundary_quart_dummy,
                                   lcrime,ses_blv,lpop,p_race_white,age_15_35_male,
                                   hhi,lmhhi,pown,p_poverty,p_emp_unemployed,pcol,BG_CODE,pop,
                                   phys_bound, com_density) %>% na.omit()

atl1 = lm_robust(lfelonies_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                   hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound +
                   com_density, 
                 data = mdata, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

summary(mdata$lcrime)
low_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(0,length(lcrime)),
  lcrime = seq(0,6.9,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

low_preds<- data.frame(cbind(prediction(atl1, low_dat, interval = "confidence", level=.95), low_dat))

head(low_preds)
hi_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(1,length(lcrime)),
  lcrime = seq(0,6.9,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

hi_preds<- data.frame(cbind(prediction(atl1, hi_dat, interval = "confidence", level=.95), hi_dat))

boundary<-rep("Boundary",length(low_preds))
ylow<-data.frame(low_preds)
ylow<-ylow %>% mutate(boundary="No Boundary")
yhi<-data.frame(hi_preds)
yhi<-yhi %>% mutate(boundary="Boundary")
preds<-data.frame(rbind(ylow,yhi))

atl_preds <- preds %>% mutate(city = "Atlanta")

## Austin

# expected values
mdata <- aus_fin %>% dplyr::select(lfelonies_sd,boundary_quart_dummy,
                                   lcrime,ses_blv,lpop,p_race_white,age_15_35_male,
                                   hhi,lmhhi,pown,p_poverty,p_emp_unemployed,pcol,BG_CODE,pop,
                                   phys_bound, com_density) %>% na.omit()


aus1 = lm_robust(lfelonies_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                   hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density, 
                 data = mdata, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

summary(mdata$lcrime)
low_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(0,length(lcrime)),
  lcrime = seq(0,8.1,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

low_preds<- data.frame(cbind(prediction(aus1, low_dat, interval = "confidence", level=.95), low_dat))

head(low_preds)
hi_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(1,length(lcrime)),
  lcrime = seq(0,8.1,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

hi_preds<- data.frame(cbind(prediction(aus1, hi_dat, interval = "confidence", level=.95), hi_dat))

boundary<-rep("Boundary",length(low_preds))
ylow<-data.frame(low_preds)
ylow<-ylow %>% mutate(boundary="No Boundary")
yhi<-data.frame(hi_preds)
yhi<-yhi %>% mutate(boundary="Boundary")
preds<-data.frame(rbind(ylow,yhi))

aus_preds <- preds %>% mutate(city = "Austin")

## Boston

# expected values #
mdata <- bos_fin %>% dplyr::select(lfelonies_sd,boundary_quart_dummy,
                                   lcrime,ses_blv,lpop,p_race_white,age_15_35_male,
                                   hhi,lmhhi,pown,p_poverty,p_emp_unemployed,pcol,BG_CODE,pop,
                                   phys_bound, com_density) %>% na.omit()

bos1 = lm_robust(lfelonies_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                   hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + 
                   com_density, 
                 data = mdata, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

summary(mdata$lcrime)
low_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(0,length(lcrime)),
  lcrime = seq(0,7,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

low_preds<- data.frame(cbind(prediction(bos1, low_dat, interval = "confidence", level=.95), low_dat))

head(low_preds)
hi_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(1,length(lcrime)),
  lcrime = seq(0,7,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

hi_preds<- data.frame(cbind(prediction(bos1, hi_dat, interval = "confidence", level=.95), hi_dat))

boundary<-rep("Boundary",length(low_preds))
ylow<-data.frame(low_preds)
ylow<-ylow %>% mutate(boundary="No Boundary")
yhi<-data.frame(hi_preds)
yhi<-yhi %>% mutate(boundary="Boundary")
preds<-data.frame(rbind(ylow,yhi))

bos_preds <- preds %>% mutate(city = "Boston")

## Chicago

# expected values
mdata <- chi_fin %>% dplyr::select(lfelonies_sd,boundary_quart_dummy,
                                   lcrime,ses_blv,lpop,p_race_white,age_15_35_male,
                                   hhi,lmhhi,pown,p_poverty,p_emp_unemployed,pcol,BG_CODE,pop,
                                   phys_bound, com_density) %>% na.omit()

chi1 = lm_robust(lfelonies_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                   hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density, 
                 data = mdata, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

summary(mdata$lcrime)
low_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(0,length(lcrime)),
  lcrime = seq(0,9.25,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

low_preds<- data.frame(cbind(prediction(chi1, low_dat, interval = "confidence", level=.95), low_dat))

head(low_preds)
hi_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(1,length(lcrime)),
  lcrime = seq(0,9.25,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

hi_preds<- data.frame(cbind(prediction(chi1, hi_dat, interval = "confidence", level=.95), hi_dat))

boundary<-rep("Boundary",length(low_preds))
ylow<-data.frame(low_preds)
ylow<-ylow %>% mutate(boundary="No Boundary")
yhi<-data.frame(hi_preds)
yhi<-yhi %>% mutate(boundary="Boundary")
preds<-data.frame(rbind(ylow,yhi))

chi_preds <- preds %>% mutate(city = "Chicago")



## Milwaukee

# expected values
mdata <- mil_fin %>% dplyr::select(lfelonies_sd,boundary_quart_dummy,
                                   lcrime,ses_blv,lpop,p_race_white,age_15_35_male,
                                   hhi,lmhhi,pown,p_poverty,p_emp_unemployed,pcol,BG_CODE,pop,
                                   com_density, phys_bound) %>% na.omit()

mil1 = lm_robust(lfelonies_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                   hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound + com_density, 
                 data = mdata, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

summary(mdata$lcrime)
low_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(0,length(lcrime)),
  lcrime = seq(0,7.5,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

low_preds<- data.frame(cbind(prediction(mil1, low_dat, interval = "confidence", level=.95), low_dat))

head(low_preds)
hi_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(1,length(lcrime)),
  lcrime = seq(0,7.5,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

hi_preds<- data.frame(cbind(prediction(mil1, hi_dat, interval = "confidence", level=.95), hi_dat))

boundary<-rep("Boundary",length(low_preds))
ylow<-data.frame(low_preds)
ylow<-ylow %>% mutate(boundary="No Boundary")
yhi<-data.frame(hi_preds)
yhi<-yhi %>% mutate(boundary="Boundary")
preds<-data.frame(rbind(ylow,yhi))

mil_preds <- preds %>% mutate(city = "Milwaukee")

## Seattle


## expected values
mdata <- sea_fin %>% dplyr::select(lfelonies_sd,boundary_quart_dummy,
                                   lcrime,ses_blv,lpop,p_race_white,age_15_35_male,
                                   hhi,lmhhi,pown,p_poverty,p_emp_unemployed,pcol,BG_CODE,pop,
                                   phys_bound, com_density) %>% na.omit()


sea1 = lm_robust(lfelonies_sd ~ boundary_quart_dummy*lcrime + ses_blv + lpop + p_race_white + age_15_35_male + 
                   hhi + lmhhi + pown + p_poverty + p_emp_unemployed + pcol + phys_bound +
                   com_density, 
                 data = mdata, cluster = BG_CODE, se_type = "stata", subset = pop > 0)

summary(mdata$lcrime)
low_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(0,length(lcrime)),
  lcrime = seq(0,7.4,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

low_preds<- data.frame(cbind(prediction(sea1, low_dat, interval = "confidence", level=.95), low_dat))

head(low_preds)
hi_dat <- mdata %>% with(expand.grid(
  boundary_quart_dummy = rep(1,length(lcrime)),
  lcrime = seq(0,7.4,.1),
  ses_blv = mean(ses_blv), 
  lpop = mean(lpop),
  p_race_white = mean(p_race_white),
  age_15_35_male = mean(age_15_35_male),
  hhi = mean(hhi),
  lmhhi = mean(lmhhi),
  pown = mean(pown),
  p_emp_unemployed = mean(p_emp_unemployed),
  pcol = mean(pcol),
  p_poverty = mean(p_poverty),
  phys_bound = mean(phys_bound),
  com_density = mean(com_density)))

hi_preds<- data.frame(cbind(prediction(sea1, hi_dat, interval = "confidence", level=.95), hi_dat))

boundary<-rep("Boundary",length(low_preds))
ylow<-data.frame(low_preds)
ylow<-ylow %>% mutate(boundary="No Boundary")
yhi<-data.frame(hi_preds)
yhi<-yhi %>% mutate(boundary="Boundary")
preds<-data.frame(rbind(ylow,yhi))

sea_preds <- preds %>% mutate(city = "Seattle")

#### stack 'em together and graph ####


## Graph
atl_preds <- atl_preds %>% dplyr::select(lcrime,fitted.fit,
                                         fitted.lwr,fitted.upr,boundary,city)

aus_preds <- aus_preds %>% dplyr::select(lcrime,fitted.fit,
                                         fitted.lwr,fitted.upr,boundary,city)

bos_preds <- bos_preds %>% dplyr::select(lcrime,fitted.fit,
                                         fitted.lwr,fitted.upr,boundary,city)

chi_preds <- chi_preds %>% dplyr::select(lcrime,fitted.fit,
                                         fitted.lwr,fitted.upr,boundary,city)

mil_preds <- mil_preds %>% dplyr::select(lcrime,fitted.fit,
                                         fitted.lwr,fitted.upr,boundary,city)

sea_preds <- sea_preds %>% dplyr::select(lcrime,fitted.fit,
                                         fitted.lwr,fitted.upr,boundary,city)

d <- rbind(atl_preds, aus_preds, bos_preds, chi_preds,
           mil_preds,sea_preds)
head(d)

p <- ggplot(data=d, aes(x=lcrime, y=fitted.fit,group=boundary)) +
  facet_wrap(~ city, nrow = 2, scales = "free")+
  #geom_point(aes(shape = Boundary), color = "grey",alpha=.5)+
  geom_line(aes(y=fitted.fit, linetype=boundary, color = boundary))+
  geom_ribbon(aes(ymin=fitted.lwr, ymax=fitted.upr,linetype=boundary, color=boundary),alpha=0.1)+
  labs(y="Felony Arrests (Logged)",x="Crime (Logged)",
       title="Police Attention to Crime in and out of Boundary Zones")+
  theme_bw(base_size=12,base_family="Times")+
  theme(panel.border=element_rect(fill=NA, colour=NA), 
        legend.title = element_blank(),
        legend.position="bottom",
        panel.grid.minor=element_line(colour=NA))+
  theme(strip.background = element_rect(fill="white"))


ggsave(plot = p, filename = "figure3_plots.png",
       height = 8, width = 8)




